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Reynolds-averaged and hybrid Reynolds-averaged/large-eddy simulations 
have been applied to a supersonic coaxial jet flow experiment. The exper- 
iment was designed to study compressible mixing flow phenomenon un- 
der conditions that are representative of those encountered in scramjet 
combustors. The experiment utilized either helium or argon as the in- 
ner jet nozzle fluid, and the outer jet nozzle fluid consisted of laboratory 
air. The inner and outer nozzles were designed and operated to produce 
nearly pressure-matched Mach 1.8 flow conditions at the jet exit. The pur- 
pose of the computational effort was to assess the state-of-the-art for each 
modeling approach, and to use the hybrid Reynolds-averaged/large-eddy 
simulations to gather insight into the deficiencies of the Reynolds-averaged 
closure models. The Reynolds-averaged simulations displayed a strong sen- 
sitivity to choice of turbulent Schmidt number. The initial value chosen for 
this parameter resulted in an over-prediction of the mixing layer spreading 
rate for the helium case, but the opposite trend was observed when ar- 
gon was used as the injectant. A larger turbulent Schmidt number greatly 
improved the comparison of the results with measurements for the helium 
simulations, but variations in the Schmidt number did not improve the 
argon comparisons. The hybrid Reynolds-averaged/large-eddy simulations 
also over-predicted the mixing layer spreading rate for the helium case, 
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while under-predicting the rate of mixing when argon was used as the in- 
jectant. The primary reason conjectured for the discrepancy between the 
hybrid simulation results and the measurements centered around issues re- 
lated to the transition from a Reynolds-averaged state to one with resolved 
turbulent content. Improvements to the inflow conditions were suggested 
as a remedy to this dilemma. Second-order turbulence statistics were also 
compared to their modeled Reynolds-averaged counterparts to evaluate the 
effectiveness of common turbulence closure assumptions. 
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Subscripts 
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inj 
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computational coordinate indices 
center jet injectant index 
coflow air index 


Introduction 

Reynolds-averaged Computational Fluid Dynamics (CFD) models have become an inte- 
gral part of the design and analysis of high-speed air-breathing engines. The maturation of 
multi-purpose CFD codes coupled with advancements in computer architectures have sub- 
stantially reduced the turn-around time required to perform steady-state Reynolds- Averaged 
Simulations (RAS). Unfortunately, the turbulence models required to close the Reynolds- 
averaged equation set have not kept pace with these advancements. As a result, RAS ap- 
proaches often require calibrations for key model parameters. One example is given in Ref. 1, 
where variations of the turbulent Prandtl and Schmidt numbers for a scramjet combustor 
simulation were shown to produce outcomes that ranged from engine unstart to complete 
flame blow-out. Similar examples that illustrate solution sensitivities to unknown RAS mod- 
eling parameters can be found in Refs. 2 and 3. 

Large Eddy Simulation (LES) methods have the potential to reduce the modeling sen- 
sitivity inherent to RAS approaches, since the intent of LES is to resolve the large scale 
turbulent structures while modeling only the small scales. Regrettably, the computational 
expense of wall-resolved LES (particularly when applied to configurations of interest to the 
high-speed propulsion community) is well beyond what can be deemed as practical by to- 
day’s standards. Hybrid RAS/LES approaches offer some relief to the computational costs 
associated with LES. These methodologies allow LES content to be resolved in areas that 
require a more rigorous modeling approach, while maintaining a more cost effective RAS 
approach for benign regions of the flow ( e.g . attached boundary layers). Hybrid approaches 
first appeared in the literature a decade ago, 4 ’ 5 and while many advancements have been 
made since then, the modeling is far from mature as evidenced by the vast array of hybrid 
approaches that have appeared in recent years. 6-9 The computational expense required for 
a hybrid simulation, while less than that of a full LES, is still formidable when compared 
with steady-state RAS. 

The present effort utilizes both RAS and hybrid RAS/LES approaches to model a series 
of supersonic coaxial jet experiments 10, 11 that have been studied at the NASA Langley 
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Research Center. These experiments involve the coaxial injection of either helium or argon 
into air. CFD simulations (based on RAS approaches) have previously been performed 
to augment the measurements in the works referenced above. Additional RAS data was 
gathered in this effort to ensure that consistent practices were employed for the comparisons 
made with hybrid RAS/LES results. The specific goals of this computational effort are to: 

• Assess the capabilities of RAS and hybrid RAS/LES for predicting high speed mixing 
layers 

• Determine the model sensitivities for both RAS and hybrid RAS/LES approaches 

• Attempt to use the hybrid RAS/LES results to assess the appropriateness of models 
used by RAS 

Measured values of composition, Pitot pressure, velocity, and rms of the velocity fluctuation 
are available for comparison with the simulations. Additional higher-order correlations (that 
are difficult to obtain experimentally) are extracted from the hybrid RAS/LES results and 
compared with the RAS predictions. The correlations considered include the Reynolds stress 
tensor and the Reynolds mass flux vector. 

Geometry Description and Flow Conditions 

A schematic of the coaxial nozzle assembly is shown in Fig. 1. The injectant supplied to 
the center jet nozzle was either a mixture of 95% helium and 5% oxygen (by volume), or pure 
argon. In the former case, a trace amount of oxygen was added to allow for the measurement 
of the streamwise component of velocity using the RELIEF 12 oxygen flow-tagging technique. 
The internal diameter of the center jet nozzle is 10 mm at the nozzle exit. The center-body 
that forms the internal nozzle is 0.25 mm thick at the nozzle exit, providing a small blunt 
base to anchor the shear layer formed between the two nozzle streams. The coflow air nozzle 
has an internal diameter of 60.47 mm, and the outer surface of this nozzle extends 12.66 mm 
downstream of the center-body. The 38.6° juncture between the internal surface of the coflow 
nozzle and the conical exterior surface is sharp. Hence, there is no appreciable base region 
to segregate the outer jet flow from the surrounding ambient air. Further details concerning 
the geometry of the rig, and the methodology used for its design, can be found in Ref. 10. 

The nominal flow conditions for each experiment are given in Tables 1 and 2 along with 
the cited uncertainties. 10, 11 The thermocouples used to measure the jet temperatures were 
located in the gas supply lines, and the pressure taps used to measure the operating pressures 
were positioned as shown in Fig. 1. Both nozzle streams have a design Mach number of 1.8. 
The flow velocity exiting the inner jet nozzle, however, is markedly different for each case, 
due to the disparate molecular weights of each center jet injectant. The velocity of the center 
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jet is more than twice that of the coflow jet for the helium case (Case 1), but it is 16 % 
lower than that of the coflow jet for the argon case (Case 2). An estimation of the convective 
Mach number: 13 

Kir ~ a) 

®air ”1“ U/ny 

yields a value of 0.7 for Case 1 and 0.16 for Case 2, respectively. Hence, one expects 
compressibility effects to be prevalent in Case 1, while Case 2 is expected to behave more 
like an incompressible shear layer. 

Computational Methodology 

All computational results were obtained using the VULCAN (Viscous Upwind aLgo- 
rithm for Complex flow ANalysis) software package. The code solves the unsteady, conser- 
vation equations appropriate for caloric-ally or thermally perfect gases with a cell-centered 
finite volume scheme. Efficient utilization of parallel architectures is realized through calls 
to MPI (Message Passing Interface) routines using an SPMD (Single Program Multiple 
Data) paradigm; a natural choice for multi-block flow solvers. Arbitrary block-to-block 
(non-aligned) connectivity is also available, allowing the flexibility to add or remove grid 
points at zonal interfaces. A variety of upwind formulations are available for evaluating the 
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Table 1: Case 1 - Helium- Air Test Conditions 


Nominal Conditions 

Center Jet 

Coflow Jet 

Ambient 

Mach Number 

1.8“ 

1.8“ 

0.025 6 

Total Temperature [K] 

305.0 (±9) 

300.0 (±6) 

294.6 (±6) 

Total Pressure [kPa] 

614.93 (±6) 

579.80 (±4) 

101.325 (±1) 


a nozzle design Mach number, b value assumed for the entrained ambient flow 


Table 2: Case 2 - Argon- Air Test Conditions 


Nominal Conditions 

Center Jet 

Coflow Jet 

Ambient 

Mach Number 

1.8“ 

1.8“ 

0.025 6 

Total Temperature [K] 

297.9 (±3.5) 

294.3 (±3.5) 

294.6 (±3.5) 

Total Pressure [kPa] 

615.86 (±5.5) 

580.68 (±4.4) 

101.325 (±0.6) 


a nozzle design Mach number, b value assumed for the entrained ambient flow 


inviscid fluxes, while central differences are used for the viscous fluxes. Diagonalized Approx- 
imate Factorization (DAF) or planar relaxation with Incomplete LU (ILU) are the primary 
options available for steady-state simulations. An explicit multi-stage Runge-Kutta or an 
implicit dual time-stepping strategy (utilizing DAF or ILU) are the options available for 
unsteady applications. A variety of one-equation and two-equation turbulence models exist 
for Reynolds-averaged simulations. LES and hybrid RAS/LES models are also available for 
turbulent flows that require a more rigorous modeling effort. Further details describing the 
code are found elsewhere. 14, 15 

RAS Numerical Model Description 

The RAS results utilized the Low-Diffusion Flux Split Scheme (LDFSS) of Edwards 16 
to evaluate the inviscid fluxes. The Monotone Upstream-centered Scheme for Conservation 
Laws (MUSCL) extrapolation parameter (/c) was chosen as 1/3 to minimize spatial trun- 
cation error, and the van Leer flux limiter 17 was employed to enforce the Total Variation 
Diminishing (TVD) property. All of the steady-state RAS solutions were advanced in time 
with a diagonalized approximate factorization scheme. 

RAS Physical Model Description 

The turbulence model chosen was the Wilcox (1998) k- u model, 18 and the wall function 
procedure of Wilcox 19 was used to relax the grid spacing requirements near solid surfaces. 
A dilatation-dissipation modification 18 to the turbulence kinetic energy equation was also 
enabled to reduce the shear layer spreading rate that has been observed under high com- 
pressible Mach number conditions. The baseline values for the turbulent Prandtl (Pr t ) and 
Schmidt ( Set ) numbers, which control the turbulent transport of energy and mass, were 
chosen as 0.9 and 0.5, respectively. The turbulent Schmidt number was one of the properties 
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that was varied in parametric studies to evaluate the sensitivity of the RAS solutions to the 
value assumed for this coefficient. 

Hybrid RAS/LES Numerical Model Description 

The hybrid RAS/LES results required a somewhat different algorithmic approach to en- 
courage the development of resolved turbulent content. In particular, the standard Total 
Variation Diminishing (TVD) limiters that are typically employed for RAS are overly dis- 
sipative when used for LES. This class of limiter tends to locally reduce the accuracy of 
the inviscid flux scheme to first-order whenever extrema of any magnitude are encountered. 
Essentially Non-Oscillatory (ENO) limiters offer a simple means of maintaining second-order 
accuracy near local extrema. This is accomplished by utilizing data at additional grid nodes 
to alter the stencil used to construct the inviscid flux terms. The ENO limiter used in this 
effort is the SONIC-A limiter of Suresh and Huynh, 20 , which can be considered as a second- 
order extension of the van Leer TVD limiter. All of the time-accurate hybrid RAS/LES 
solutions were advanced in time using a dual time-stepping approach that combined the di- 
agonalized approximate factorization scheme for integration in pseudo-time, with a 3-point 
backwards finite difference approximation for integration in real-time. The values selected for 
the physical time-step and sub-iteration CFL constraint were 0.05 [is and 50.0, respectively. 
The time-step was chosen based on cell residence time considerations to ensure that turbu- 
lent structures would not traverse more than one grid cell per time-step. The sub-iteration 
process was considered converged when the residual error dropped at least 2.5 orders of 
magnitude. This level of convergence typically required 4-7 sub-iterations for each physical 
time-step. 

Hybrid RAS/LES Physical Model Description 

The hybrid RAS/LES methodology used in this effort builds on the previous work discussed 
in Ref. 6. This framework was designed to enforce a RAS behavior near solid surfaces, and 
switch to an LES behavior in the outer portion of the boundary layer and free shear regions. 
Hence, this formulation can be thought of as a wall-modeled LES approach, where RAS is 
used as the near-wall model. The basic idea is to blend any trusted RAS eddy viscosity with a 
desired LES Sub-Grid Scale (SGS) viscosity, along with any transport equations that involve 
a common RAS and SGS property. In this effort, the Wilcox (1998) k-u RAS model 18 was 
blended with the one-equation SGS model of Yoshizawa. 21 The Yoshizawa model involves 
an evolution equation for the SGS turbulence kinetic energy, hence the blended expressions 
that are appropriate for this model combination are: 

Hybrid RAS/SGS viscosity = (F) [RAS viscosity] + (1 — F) [SGS viscosity] 

Hybrid RAS/SGS k-equation = (F) [RAS k-equation] + (1 — F) [SGS k-equation] (2) 
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where F is a blending function that varies between 0 and 1. Note that the transport equation 
for the RAS specific dissipation rate (u) does not have an SGS counterpart. Hence, the 
blending is not applied to this equation, and all of the terms in this equation that involve 
the eddy viscosity are evaluated based on RAS relationships. 

The motivation behind the development of this hybrid RAS/LES framework is two-fold. 
First, the blending of two independent RAS and LES closure models offers the flexibility of 
having an optimized set of closure equations for both RAS and LES modes. The second (and 
more critical) driving factor was the desire to alleviate the difficulties associated with the 
design of grid topologies that are appropriate for purely grid-dependent blending paradigms 
such as those used for Detached Eddy Simulation (DES). 4i 22 The movement away from 
simple grid-dependent blending strategies has gained momentum in recent years. 7) 8 In 
fact, even the developers of DES are now promoting an “improved” version of their scheme 
termed Delayed Detached Eddy Simulation (DDES) 8 that involves flow-dependent functions 
to switch between RAS and LES regimes. 

The blending function (F) used in this effort is based on the ratio of the wall distance d 
to a modeled form of the Taylor microscale (£): 



where k is the von Karman constant (0.41), C M is 0.09, a is a user-defined model constant, 
and (j) is set to tanh~ l (0.98) to force the balancing position of F (be. the position where 
nrj 2 = \J~C~p) to 0.99. The value chosen for a provides control over the y + position where the 
average LES to RAS transition point (defined as F — 0.99) occurs. If resolved LES content 
is desired for an attached boundary layer, then this constant should be set such that the 
transition point occurs in the region where the boundary layer wake law starts to deviate 
from the log law. If the transition point is enforced at a lower y + value (e.g. well within 
the log law region), a dual log layer can appear. 23, 24 Conversely, if the transition point 
is enforced at a y + value that extends well into the defect layer, then the level of resolved 
turbulence will be limited. Details on a procedure to analytically determine the value for a 
that corresponds to a target y + value is described in Ref. 24. 
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Grid and Boundary Condition Details 

A two-dimensional (axisymmetric) grid was generated for the steady-state Reynolds- 
averaged solutions (see Fig. 2). This grid consisted of just under 250,000 cells divided across 
five structured grid zones. The computational domain extended 150 center jet diameters 
(D) downstream of the jet exit, and between 40 and 70 jet diameters in the radial direction. 
The internal flowfield within the concentric nozzles was also included as part of the com- 
putational domain to allow for the development of boundary layers. The last experimental 
data stations were located 26.101 (Case 1) and 45.276 (Case 2) jet diameters downstream of 
the center jet exit plane. Hence, most of the grid was clustered between the jet exit and an 
x/D of 50. The portion of the computational domain downstream of x/D = 50 was added to 
provide plenty of space for the jet flow to reach a subsonic state, which ensured that the spec- 
ified pressure outflow boundary condition (with all other variables extrapolated) remained 
well-posed. The far-held boundary was also placed many tens of jet diameters away from 
the domain of interest to minimize any chance of data corruption that could occur via wave 
reflections off of the characteristic inflow condition 25 that was applied along this boundary. 
The nozzle stagnation conditions given in Tables 1 and 2 were applied at each nozzle inflow 
plane, with the Mach number extrapolated from the interior to complete the specification of 
these subsonic inflow boundaries. The operating total temperature of each nozzle was nom- 
inally the same as the room temperature, so all solid surfaces were assumed to be adiabatic 
no-slip surfaces. The grid was clustered to all solid surfaces at a level appropriate for the 
use of wall functions ( y + < 36). Finally, the entire grid system was split up into a total of 
179 grid blocks, which yielded excellent load balance statistics for up to 64 processors. 

An azimuthal slice of the three-dimensional grid generated for the hybrid RAS/LES cases 
is shown in Fig. 3. The boundary conditions and the extent ( x and r) of the computational 
domain was identical to that used for the axisymmetric RAS cases. The gridding strategy, 
however, was altered to reflect the change in computational algorithm. LES requires the 
use of grid cells that are roughly isotropic to allow for the development and sustainment of 
resolved turbulent structures. Hence, the streamwise grid spacing had to be reduced from 
that utilized for the pure RAS. This requirement, coupled with the fact that the full 3-D 
flowfield must be solved, forced a concession to be made to keep the computational costs 
within reason. The compromise that was accepted was to provide a level of grid resolution 
capable of supporting LES only up to an x/D of 25 (highlighted in red). This domain 
captures all but the last experimental station for Case 1, and 13 of the 16 stations for 
Case 2. Non-aligned zonal interfaces (patches) were inserted at x/D stations of 25, 50, 
and 100 to systematically remove grid nodes in the radial direction. A patch interface was 
also inserted in the ambient air region to coarsen the streamwise spacing in the far-held 
where axial refinement is not required. Each patch interface is denoted by a solid black 
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Figure 2: Axisymmetric grid utilized for RAS (247368 cells) 


§ 



Figure 3: Azimuthal slice of the three-dimensional grid utilized for hybrid RAS/LES 

(43,285,632 cells) 



Figure 4: 


Isometric visualization of the LES-resolved portion of the hybrid RAS/LES grid 
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line in Fig. 3. In addition to substantially reducing the computational costs, the coarsening 
process performed at large x/D tends to dissipate the large scale vortical structures that 
propagate out of the LES domain of interest. This reduces the likelihood of an ill-posed 
boundary condition appearing at the outflow plane. An H-0 mesh topology was chosen 
for the cross-flow planes to avoid the numerical difficulties associated with collapsed faces 
that would result from a classic polar topology. This feature of the grid is illustrated in the 
three-dimensional view of the LES domain shown in Fig. 4. Finally, the entire grid system 
(43,285,632 cells) was split up into a total of 1669 grid blocks, which resulted in excellent 
load balance statistics for up to 360 processors. 

Reynolds- Averaged Results 

Contours of helium mass fraction and Mach number extracted from the baseline RAS 
(, Sc t = 0.5) are shown in Fig. 5. The helium mass fractions have been normalized by the the 
center jet value (95% by volume corresponds to 70.39% by mass): 


Y He => WO. 7039 (4) 

for all of the results that follow to simplify the analysis and comparisons with the argon 
cases. The helium contours show that the potential core of the center jet persists for ap- 
proximately 12 jet diameters, and the center-line helium mass fraction in the plume that 
forms downstream of this station decays to a value of 0.115 at an x/D = 50. The Mach 
contours show that both nozzle streams are nearly pressure matched. Although not visible at 
this scale, a pair of counter-rotating separation bubbles are present behind the small blunt 
base of the center-body. The weak expansion fan that is seen on either side of the base 
region forms as the jets negotiate around this separated flow zone. A weak recompression 
shock then develops as the flow is forced to turn back on itself at the close-off point of the 
recirculation bubble. The conical recompression shock that formed in the center jet flowheld 
steepens as it approaches the axis of symmetry, resulting in a Mach disk at the axis. As will 
be shown later, the total pressure loss that occurs across this shock wave is also evident in 
the measured Pitot pressure surveys. 

Measured helium mass fraction profiles are compared with the RAS results in Fig. 6. The 
measurements of injectant composition were made using a gas sampling probe as detailed 
in Ref. 10. In general, the Set = 0.5 simulation over-predicted the rate of mixing between 
the coaxial streams. Hence, an additional simulation was performed with the Sc t value 
increased by a factor of 2 (to 1.0). The third result shown in this figure was obtained from 
a coarse grid simulation (factor of 2 coarsening in each coordinate direction) that used an 
Sc t value of 0.5. Comparisons are made at 8 axial stations downstream of the center jet 
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Y He : 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.50 0.55 0.70 0.75 0.80 0.B5 0.90 0.95 1.00 



Mach: 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.50 1.70 1.B0 1.90 2.00 


Figure 5: Contours of normalized helium mass fraction and Mach number for the baseline 

RAS 

exit plane. The results extracted from the Sc t = 0.5 simulation over-predicted the measured 
growth rate of the mixing layer for all but the first experimental station. This feature 
is illustrated by the wider mixing width between the two streams (as compared with the 
measurements), and leads to a premature close-off of the potential core of the helium jet. 
It should be noted that the RAS model used in this effort included a dilatation-dissipation 
modification 18 to reduce the shear layer spreading rate for compressible mixing layers. This 
leaves the turbulent Schmidt number (which controls the rate of turbulent mass diffusion) 
as the remaining parameter that can be adjusted to reduce the mixing rate. The results 
obtained when doubling the Sc t value (which reduces the modeled turbulent diffusion term 
by a factor of 2) agree well with the measurements. In fact, the results obtained with 
Set = 1.0 matched nearly every measured result to within the precision of the measurement 
device of ±1.0 — 1.5%. The effect of coarsening the grid resolution had a much smaller 
impact on the predicted results, indicating that the level of grid convergence is well within 
the uncertainty bounds of the physical models employed. 

Measured Pitot pressure profiles 10 are compared with the RAS results in Fig. 7. The 
large deficit in Pitot pressure near r/D = 0.5 is a result of viscous effects from the upstream 
boundary layers. The small deficit seen near the axis of symmetry at the first four axial 
stations is a result of the Mach disk formed by the coalescence of the conical shock front 
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that was discussed previously (see Fig. 5). The measured value of this deficit is slightly 
larger than the predicted value. The Pitot pressure deficit is also under-predicted for the 
coflow shock structure that appears near r/D = 1.1 at the first axial station. In general, the 
best agreement with measurements is given by the Sc t = 1.0 data. The level of agreement 
seen between this simulation and the measured values is within the measurement uncertainty 
(based solely on the transducer error) of ±0.5% at nearly every measured point for x/D < 15. 
Downstream of this station, the Sc t = 1.0 Pitot data under-predicts the measured values 
near the axis, and over-predicts them near the lower edge of the shear layer present between 
the coflow jet and the ambient surroundings. Once again, the differences noted between the 
fine and coarse grid solutions is much smaller than that associated with a modification of 
the turbulent Schmidt number. 

Contours of argon mass fraction and Mach number extracted from the RAS (Sc t = 0.5) 
are shown in Fig. 8. The argon results show a much longer potential core than that for the 
helium case. In fact, the potential core in the simulations extends beyond the x/D = 50 
station shown in Fig. 8. The slower mixing rate between the two jet streams is a result of 
the reduced shear associated with Case 2. The velocity of the argon stream is only 16% 
lower than the coflow stream. This is contrasted with the helium case, where the velocity 
of the helium stream was more than twice that of the coflow stream. The Mach contours 
clearly show that the nozzle streams are nearly pressure matched, and the near-field flow 
structures are practically identical to those discussed previously for Case 1. The far-held 
shock/expansion pattern is quite different, however, due to the reduced spreading rate of the 
mixing layer. 

Measured argon mass fraction profiles 11 are compared with the RAS results in Fig. 9. The 
Sc t = 0.5 argon simulation under-predicted the rate of mixing between the coaxial streams. 
This is in contrast to the helium case where the same value of turbulent Schmidt number 
over-predicted the growth rate of the mixing layer. Hence, an additional simulation was 
performed with the turbulent Schmidt number reduced by a factor of 2. A third simulation 
was performed (with Sc t = 0.5) using the coarse grid to evaluate the grid dependence. 
Comparisons are made at the same 8 axial stations used to evaluate the helium injection 
data. The results obtained when reducing the Sc t value had the desired effect of reducing 
the length of the potential core of the argon jet. However, this adjustment also resulted in a 
peculiar inflection point near the outer edge of the mixing layer that was not captured in the 
measurements. The reduced Schmidt number results also over-predicted the spreading rate of 
the mixing layer for x/D < 18. Overall, the best agreement with measurements was obtained 
with an Sc t value of 0.5, which predicted the correct spreading rate of the shear layer for 
x/D < 18. The convective Mach number for the mixing layer was estimated to be 0.16, so an 
additional simulation (not shown) was performed with the dilatation-dissipation modification 
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* Measured 

RA5|Sc,=0.5) 

RA5 (Sc,= 1.0) 

RA5 (Sc-0.5, coarse) 










Figure 6: Comparison of normalized helium mass fraction (RAS predictions) with measured 
values 
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* Measured 

RAS (Sc, =0.5) 

RAS (Sc,= 1.0) 

RA5 (Sc -0.5, coarse) 










Figure 7: Comparison of Case 1 Pitot pressure (RAS predictions) with measured values 
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Y Ar : 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.B5 0.90 0.95 1.00 



Mach: 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.B0 1.90 2.00 

Figure 8: Contours of argon mass fraction and Mach number for the baseline RAS 

disabled to ensure that the model was not impacting the predictions. As expected, the 
compressibility correction had no impact on this simulation. One interesting feature noticed 
in the experimental data is an increased deviation from symmetry that begins to appear 
at the x/D = 12 station. This happens to be the station where the simulations begin to 
deviate from the measurements. The experimental probes were traversed over the entire 
coaxial flowheld (i.e. data was gathered for both positive and negative ^/-values relative to 
the jet axis). Since the data should be axisymmetric in the mean, the measured data was 
plotted against the absolute value of y in this work, which effectively results in two sets of 
“radial” data at each experimental station. The deviation from symmetry is much larger for 
Case 2 than for Case 1, which might have implied a slight misalignment of the jet assembly 
during the argon tests. 11 As was seen for the helium injection conditions, the predictions 
offered by both the fine and coarse meshes were quite similar. 

Measured Pitot pressure profiles 11 for Case 2 are compared with the RAS results in 
Fig. 10. The near-held Pitot profile features are identical to those seen for Case 1, as implied 
by the Mach contours that were discussed earlier. In general, the level of agreement seen 
between the simulations and the measured values is within the measurement uncertainty for 
x/D < 15. Downstream of this station, the computed Pitot data shows sharper features, 
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which is a direct result of the under-prediction of the mixing processes at this point in the 
flowfield. It is interesting to note that the over-prediction of the Pitot pressure near the 
lower edge of the shear layer present between the coflow jet and the ambient surroundings 
(r/D ~ 2.25) is larger that what was seen for Case 1. The shear between the coflow jet and 
the ambient environment is practically identical to that for the helium condition. Hence, the 
under- prediction of the argon/air mixing layer has evidently modified the shock system in 
the outer jet enough to alter the behavior of the coflow/ ambient air shear layer. 


Hybrid RAS/LES Results 


The coaxial flowfields examined in this effort have three distinct flow regions: center jet, 
coflow jet, and the ambient surroundings. The core velocity in the coflow jet was 480 m/s, 
and the center jet values were 1100 m/s (helium) and 400 m/s (argon). The ambient 
air region was almost stagnant. This presented a problem in that the hybrid RAS/LES 
approach requires a time-accurate integration with a time-step small enough to temporally 
resolve turbulent structures at scales representative of a computational cell width. A time- 
step appropriate for resolving a turbulent structure in the jet flow is given by the following 
estimate: 


At 


Ax 


TflClX [Ui n j , Hair ) 
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Based on the disparate characteristic velocities between the jet flows and the ambient air, 
this definition would require hundreds of time-steps to transport a particle across one cell 
length in the ambient air region. Fortunately, the shear layer of interest is the mixing layer 
between the supersonic jets rather than the shear layer present between the coflow jet and the 
ambient air. Hence, the time scale dilemma was dealt with by forcing the hybrid RAS/LES 
model to remain in RAS mode for the outer portion of the coflow jet and ambient region. In 
particular, an initial solution was obtained for the entire flowfield using a steady-state RAS. 
The calculation was then restarted using the hybrid RAS/LES approach with the blending 
function forced to unity for r/D > 2. This value was found to be large enough to contain 
all convecting coherent structures formed in the mixing layer between the two jet flows for 
the LES domain of interest. 

The form of SGS closure used in this effort was based on the one-equation formulation of 
Yoshizawa. 21 The Yoshizawa model involves an evolution equation for the SGS turbulence 
kinetic energy: 
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Figure 9: Comparison of argon mass fraction (RAS predictions) with measured values 
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Figure 10: Comparison of Case 2 Pitot pressure (RAS predictions) with measured values 
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However, two noteworthy modifications to the model were made for this effort. The first 
modification involves the definition of the filter width, A. Traditional LES practitioners have 
tended to work with truly isotropic grids, justifying the use of a simple isotropic definition 
for the filter width ( e.g . the cubed root of the cell volume). Grids utilized for configurations 
of engineering interest, however, will have some level of anisotropy. Hence, hybrid RAS/LES 
practitioners tend to prefer the following anisotropic definition for the filter width: 


A = max (Aj, A j, A*,) (8) 

In this expression, Aj, Aj, and A& denote the width along each computational direction 
of a three-dimensional hexahedral grid cell. The second modification affects the level of 
SGS viscosity in LES regions. At equilibrium (defined when production balances dissipation 
in Eq. 6), the Yoshizawa model reduces to an algebraic Smagorinsky closure model. The 
effective Smagorinsky constant (C s ) implied by the model under equilibrium conditions is 
given by the following expression: 

(C s f=V2C l ,( < Py (9) 

The default values given in Ref. 21 for C and C r j are 0.05 and 1.0, respectively. The value 
of C s implied by these values is 0.126. There is no universal value for this Smagorinsky 
constant, but the LES community has typically endorsed values on the order of 0.2 for 
homogeneous turbulence and 0.065 for shear flows. 26 Initial calculations used a reduced value 
of 0.02075 for C' /t to recover the “accepted” Smagorinsky constant for shear flows. This value, 
however, resulted in a solution with an excessive delay in the formation of resolved turbulent 
structures. The numerical dissipation of the second-order upwind numerical framework used 
in this effort is on the same order as the SGS viscosity (A 2 ), hence one might expect lower 
levels of SGS viscosity to be required. A value of = 0.00823 was found to readily 
unlock significant levels of resolved turbulence energy in the present simulations. This value 
enforces a Smagorinsky constant of 0.0325 (half the recommended value for shear flows) 
under equilibrium conditions, and was selected as the baseline setting. 

An instantaneous snapshot of the flowfield for Case 1 is shown in Fig. 11. The top im- 
age displays the density gradient magnitude (numerical approximation of a Schlicren image) 
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Figure 11: Instantaneous Schlieren and normalized helium mass fraction contours (hybrid 

RAS/LES) 
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to illustrate the flow structure, and the bottom image shows the normalized helium mass 
fraction contours. A significant level of resolved turbulent content was captured in the sim- 
ulation. The recirculation zone behind the base of the center-body provided a source of 
large scale unsteadiness that rapidly triggered Kelvin Helmholtz instabilities in the shear 
layer. These instabilities transitioned the flow from a fully modeled Reynolds-averaged state 
exiting the nozzle, to a primarily resolved turbulent state a few jet diameters downstream of 
the nozzle exit plane. Unsteady shock and expansion waves are seen reflecting through the 
jet structure. The shear layer between the outer jet and ambient air maintained the desired 
steady-state Reynolds-averaged behavior due to the forced usage of the RAS equations out- 
side of r/D = 2.0. The high level of shear associated with this operating condition provided 
a strong mechanism for turbulence production as illustrated by the rich array of vortical 
structures seen in the helium contours. These structures provide the turbulent “stirring” 
that leads to enhanced mixing by stretching the interface between the air and injectant. 

Prior to gathering any statistics, time integrations were performed for two complete 
flow-through times to flush out the initial transients during the transition from the con- 
verged RAS solution to the hybrid RAS/LES result. This relatively small time interval was 
sufficient due to a lack of any large regions of separated flow, and the fact that the am- 
bient region was essentially frozen in its Reynolds-averaged state. At this point, running 
averages were computed and stored after each time-step. In order to reduce the integration 
time required to gather meaningful statistics, a combination of temporal and spatial (in the 
circumferential direction) averaging was performed. The spatial-average was made possible 
by a three-step process. The first step was to interpolate the ensemble-averaged data onto a 
polar grid topology (in the cross-flow plane) at the experimental streamwise stations. The 
second step involved a coordinate transformation of the flow data to a cylindrical coordinate 
system (x,y,z — > x,r,9). The final step was to average the transformed variables along 
the azimuthal curves of constant radius (which corresponds to the spatially homogeneous 
direction). All three steps were performed using the TECPLOT post-processing package 27 
and automated using the TECPLOT macro scripting language. 

The averaged hybrid RAS/LES helium mass fraction profiles are compared with the 
measurements in Fig. 12. In general, the predictions show a shear layer growth rate that 
is more rapid than that implied by the measurements. As mentioned previously, the SGS 
model coefficients were adjusted to promote the onset of flow instabilities. It appears that 
this adjustment leads to an inappropriate level of modeled viscosity (be. a value that is too 
low) once the Kelvin Helmholtz instabilities have transitioned to a fully turbulent state. This 
situation might be improved if resolved turbulent structures were present at the nozzle exit 
plane. This would alleviate the delay in transition that is associated with the conversion of 
modeled Reynolds-averaged turbulence to resolved LES content, and would allow a larger 
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Smagorinsky coefficient to be utilized. One reasonably economic approach to this task is to 
implement recycling/ rescaling techniques. 28, 29, 30 This will be the focus of future efforts. A 
comparison of the averaged hybrid RAS/LES results with those obtained using pure RAS 
at the default (Sc t = 0.5) condition show many similarities. This particular RAS simulation 
also over-predicted the shear layer growth rate, but one noteworthy difference between the 
results is the manner in which the helium mass fraction asymptotically approaches zero at the 
edge of the mixing layer. RAS approaches are notorious for producing a sharp non-physical 
interface at the edge of the shear layer. The hybrid RAS/LES results show a much smoother 
asymptotic decay that is more representative of the measurements. Statistical convergence 
has been assessed by comparing the statistics extracted over successively larger time-step 
intervals. The statistics extracted over twenty and thirty thousand time-steps (cycles) show 
practically identical results. Another useful statistical convergence check is to examine the 
averaged azimuthal velocity component. This value should approach zero as the statistics 
are converged. The maximum absolute value of azimuthal velocity (averaged over 30,000 
cycles) was 2.5 m/s, or roughly 0.2% of the peak streamwise velocity. 

The averaged Pitot pressure profiles are compared with the measurements in Fig. 13. 
These results reinforce the assertion that the mixing layer growth rate has been over- 
predicted, particularly at stations downstream of an x/D of 4. Although the deficit in 
Pitot pressure profile is broader across the shear layer, the qualitative features of the profile 
shape mimics that of the measurements. This was not the case in the RAS results that 
over-predicted the level of mixing. The finer details of the flow structure, e.g. the slight 
pressure deficit due the the Mach disk at the axis, are also evident in the hybrid RAS/LES 
data. As was the case with the helium mass fraction statistics, the averaged Pitot pressure 
profiles for both averaging intervals show practically identical results, suggesting that the 
first-order moments are well converged. 

An instantaneous snapshot of the flowheld for Case 2 is shown in Fig. 14. The top 
image displays a numerical Schlieren of the flow structure, while the lower image displays 
the argon mass fraction distribution. The Kelvin Helmholtz instabilities are clearly visible 
in both images and persist for many jet diameters before breaking down into fully three- 
dimensional turbulent structures. This transition process begins near an x/D of 10 or 12, 
where the last visible wave structures appear in the Schlieren image. The low level of shear 
for this case (as compared with Case 1) is the root cause for the delay in the break-down of 
the Kelvin Helmholtz instabilities. These instabilities are also primarily two-dimensional in 
nature; an expected result given the low convective Mach number of this mixing layer. Large 
(be. on the order of the jet diameter) turbulent coherent structures do not appear until an 
x/D of 15 to 20. Hence, the turbulent “stirring” processes that lead to enhanced mixing are 
not prevalent until the latter stages of the LES domain. 
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Figure 12: Comparison of normalized helium mass fraction (averaged hybrid RAS/LES pre- 
dictions) with measured values 
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Figure 13: Comparison of Case 1 Pitot pressure (averaged hybrid RAS/LES predictions) with 
measured values 
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Figure 14: Instantaneous Schlieren and argon mass fraction contours (hybrid RAS/LES) 


26 of 36 


The averaged hybrid RAS/LES argon mass fraction profiles are compared with measure- 
ments in Fig. 15. Overall, the level of turbulent mixing between the argon and air jets was 
under-predicted. These results are comparable with those obtained using the RAS approach 
(see Fig. 9), though the hybrid RAS/LES shows a slightly reduced mixing level at the last 
experimental station. The smooth asymptotic transition of the argon profile towards zero at 
the edge of the mixing layer is again captured in the hybrid RAS /LES results. The reduced 
spreading rate of the shear layer is also seen in the Pitot pressure profile comparisons dis- 
played in Fig. 16. The peak shear layer deficit is consistently over-predicted at each station. 
Future work will examine the influence of providing realistic resolved turbulent content at 
the exit of the coaxial nozzles. This may be particularly important for the present applica- 
tion because the boundary layer thickness of the approach coflow is an order of magnitude 
thicker than the blunt base region that is currently triggering the unsteadiness. The statis- 
tics gathered after 10,000 cycles show practically identical results as the data averaged over 
20,000 cycles, suggesting that the statistics are converged for these first-order correlations. 
The maximum azimuthal velocity (averaged over 20,000 cycles) was 0.6 m/s, or roughly 
0.1% of the peak streamwise velocity. 


Second-Order Correlations 


The primary second-order correlations that are modeled in the RAS equation set are the 
Reynolds stress tensor and the Reynolds heat and mass flux vectors. The Reynolds stress 
tensor is typically modeled by the Boussinesq approximation: 
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The accuracy of the modeling employed for these terms is critical for the success of RAS 
approaches. Measured data for each of the above correlations is scarce, so it is often difficult 
to directly assess the accuracy of the RAS closures. One of the primary goals of this effort 
was to use the hybrid RAS/LES simulations to assess the performance of the RAS modeling. 
Ideally, one would hope for a close match between the hybrid RAS/LES results and available 
measurements, but the level of agreement obtained in this effort was not entirely satisfying. 
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Figure 15: Comparison of argon mass fraction (averaged hybrid RAS/LES predictions) with 
measured values 
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Figure 16: Comparison of Case 2 Pitot pressure (averaged hybrid RAS/LES predictions) with 
measured values 
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Nevertheless, some insight into the RAS modeling can be gleaned by comparing the modeled 
terms with the values extracted from the resolved LES held. 

The rms of the streamwise velocity fluctuation is considered first. This quantity was the 
only second-order correlation that was measured, and the measurement was only taken for 
Case 1. These measurements were made using the RELIEF 12 oxygen flow-tagging technique. 
Figure 17 compares the measured values with the values computed from the RAS (Sc t = 0.5) 
result using Eq. 10, as well as the ensemble-averages extracted from the hybrid RAS/LES. 
The modeled RAS values agree remarkably well with the measurements in the early stages of 
the shear layer development, but the width of the profile is under- predicted at stations further 
downstream. The hybrid RAS/LES results over-predict the rms values in the shear layer 
at the first station, but the profiles at the stations further downstream compare favorably 
with the measurements. One particular noteworthy feature is the build-up of turbulent 
fluctuations in the core flow near the axis that is present in the hybrid RAS/LES results. 
This feature is absent from the RAS result because the mean flow gradients are relatively 
small outside of the shear layer. Mean velocity gradients are the sole source of turbulence 
production for most RAS models. The enforcement of a RAS behavior for r/D > 2 in the 
hybrid RAS/LES simulations prevented the formation of unsteady turbulence structures in 
the outer jet/ambient air shear layer; limiting the build-up of rms velocity levels in the core 
of the outer jet. Hence the core flow values in the outer jet, while larger than that predicted 
by pure RAS, was under- predicted. 

The rms radial and azimuthal velocity fluctuations are compared in Figs. 18 and 19. 
Linear eddy viscosity models based on the Boussinesq approximation tend to produce nearly 
isotropic velocity variances (normal stresses) when applied to strain dominated flows. The 
degree of anisotropy is governed by the magnitude of the velocity gradients (other than those 
contributing to the strain rate) relative to the 2/3 pk term as shown by Eq. 10. Hence, it 
is not surprising that the RAS model returned a nearly isotropic set of normal stresses for 
this turbulent shear flow. The hybrid RAS/LES data, on the other hand, showed a non- 
negligiblc level of anisotropy in the velocity variances. In particular, the streamwise velocity 
variance was substantially larger than the cross-flow variances. This result is typical for 
shear dominated flows. 

The Reynolds shear stress (u'v") and Reynolds mass flux Y^ e v" terms are compared in 
Figs. 20 and 21. The accuracy of the RAS equations, to a large extent, is driven by how 
well these terms are modeled. Since the Sc t = 1.0 RAS data compared favorably to the 
measured mean flow properties, correlations computed from this simulation were also added 
to the figures. The shear stress extracted from the hybrid RAS/LES data is larger than 
those produced by the RAS models at the first two axial stations. The situation is reversed 
at the last two stations. This result is consistent with the fact that the hybrid RAS/LES 
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results predicted a more rapid mixing rate than the RAS approaches. A similar behavior 
is evident in the Reynolds mass flux vector. The hybrid RAS/LES data showed the most 
rapid mixing, the Sc t = 1.0 RAS yielded the slowest mixing rate, and the Sc t = 0.5 RAS 
result was somewhere inbetween. Hence at the early axial stations, the turbulent mixing 
rate should be largest for the hybrid RAS/LES and smallest for the Sc t = 1.0 RAS. Once 
the mixing process has reached some critical level, the mixing rate will begin to decay. Thus 
at some station downstream, a simulation with a slower mixing rate should exceed one with 
a more rapid mixing rate. This is precisely the situation that occurs between the second and 
third streanrwise stations shown in Fig. 21. 




Figure 17: Comparison of Case 1 streamwise velocity fluctuation rms with measured values 
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Figure 18: Comparison of Case 1 radial velocity fluctuation rms 






Figure 19: Comparison of Case 1 azimuthal velocity fluctuation rms 
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Figure 20: Comparison of Case 1 u v" (Reynolds shear stress correlation) 






Figure 21: Comparison of v” (Reynolds mass flux correlation) 
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Summary 


Reynolds- averaged and hybrid Reynolds-averaged/large-eddy simulations have been per- 
formed to model a supersonic coaxial jet flow experiment. The experiment consists of an 
outer jet of air, and an inner jet that was either a He — O 2 mixture or pure argon. Both 
jets exhausted into an ambient environment. The Mach 1.8 nozzle flows exiting the test 
apparatus were nearly pressure-matched for both injectant conditions. However, the level of 
shear and the compressibility of the mixing layer varied depending on the injectants involved. 
The helium condition resulted in a highly compressible mixing layer with a convective Mach 
number of 0.7. The argon condition, on the other hand, was nearly velocity- matched and 
produced a mixing layer with a convective Mach number of 0.16. A comprehensive set of 
measurements were taken which include: Pitot pressure, mean and rms velocities, and gas 
sampling. The model geometry, flow conditions, and measurement uncertainties were all 
well documented, resulting in a package that is well suited for model validation efforts. 

The goal of the computational effort was to assess the state-of-the-art for both RAS and 
hybrid RAS/LES approaches as applied to compressible turbulent flows. The Reynolds- 
averaged simulation for the helium cases displayed a strong sensitivity to choice of turbulent 
Schmidt number. A value of 1.0 was found to be an optimal choice for this flow condition. 
A lower turbulent Schmidt number, on the order of 0.5, provided the best match with mea- 
surements for the argon case, however. The uncertainty involved with appropriate choices 
for RAS modeling parameters highlights the difficulty with using these approaches for pre- 
dictive simulations. In principle, LES or hybrid RAS/LES methods have the potential to 
reduce this uncertainty by resolving a substantial fraction of the turbulent flowfield. The 
hybrid simulations performed in this effort, however, were no more predictive than the base- 
line Reynolds-averaged predictions. The explanation provided for the discrepancy between 
the hybrid RAS/LES results and the measurements centered around issues related to how 
the Reynolds-averaged state transitions to a state with resolved turbulent content. The ad- 
dition of resolved turbulent content to the inflow conditions was suggested to address this 
concern. Finally, comparisons made between resolved second-order turbulence statistics and 
their modeled Reynolds-averaged counterparts were discussed. Refined hybrid simulations 
are required, however, to improve the accuracy of the hybrid RAS/LES results before any 
solid conclusions can be drawn concerning the RAS modeling. 
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